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ABSTRACT 



Context. Supernova remnants (SNRs) are thought to be the main source of Galactic cosmic rays (CRs) up to the "knee" in CR 
spectrum. During the evolution of a SNR, the bulk of the CRs are confined inside the SNR shell. The highest-energy particles leave 
the system continuously, while the remaining adiabatically cooled particles are released when the SNR has expanded sufficiently and 
decelerated so that the magnetic field at the shock is no longer able to confine them. Particles escaping from the parent system may 
interact with nearby molecular clouds, producing y-rays in the process via pion decay. The soft gamma-ray spectra observed for a 
number of SNRs interacting with molecular clouds, however, challenge current theories of non-linear particle acceleration that predict 
harder spectra. 

Aims. We study how the spectrum of escaped particles depends on the time-dependent acceleration history in both Type la and core- 
collapse SNRs, as well as on different assumptions about the diffusion coefficient in the vicinity of the SNR. 
Methods. We solve the CR transport equation in a test-particle approach combined with numerical simulations of SNR evolution. 
Results. We extend our method for calculating the cosmic-ray acceleration in SNRs to trace the escaped particles in a large volume 
around SNRs. We calculate the evolution of the spectra of CRs that have escaped from a SNR into a molecular cloud or dense shell 
for two diffusion models. We find a strong confinement of CRs in a close region around the SNR, and a strong dilution effect for CRs 
that were able to propagate out as far as a few SNR radii. 
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Supernova remnants (SNRs) are now widely considered to 
be sources of Galactic cosmic r ays (CRs). Diffusive shock 
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acceleration (PSA) dAxford et alJll977t iKrvmskiil Il977b iBel 
1978t lBlandford &Ostrikerlll978l) and its n on-linear modifi- 
cation (NDSA) (Malko v & O'C Drurvl |200T1) predict a power- 
law distribution (N(E) cc E~ s ) of relativistic particles with 
index s = 2 (s < 2 for the high-energy tail in NDSA). 
Despite recent advances in DSA and NDSA, it is still un- 
clear why we o bserve soft (s > 2) CR spectra at Earth 
(lAve etalJl2009h and soft ga mma-ray spectra from a number 
of SNRs: RX J0852.0- 4622 dAharonian et alj|2007b. R CW 86 
(lAharonian etaD l2009h. SN 1006 (lAcero et all |2010Y Cas A 



2010; Abdo et al.1 12010; 
l201lHGiordano etalj|201 



and Tycho's SNR 



Acciari et al 
dAcciari et al 

Understanding gamma-ray SNR spectra in various situations 
requires the study of the acceleration of the particles in the sys- 
tem (consisting of one or two shock waves), followed by the es- 
cape of these particles, their diffusion into the Galactic medium 
(Ptus kin et al.l 1201 Ot ICaprioli et al.l 120101) . and computation of 
the y-ray emission by various processes. In an earlier paper 
dTelezhinskv et al.|l2012l) . we studied the acceleration of parti- 
cles at SNR shock fronts. In this paper, we follow up this work 
by studying the escape of these particles from the system, and 
their interaction with dense media. 

A useful probe of the escape of particles from SNRs may 
be high-energy gamma-ray emission from molecular clouds lo- 
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cated close to an efficiently accelerating SNR, the so-called 
MC-SNR systems. A number of these have been observed by 
Cherenko v telescopes (CTs) and F ermi. Most sources , includ- 
ing W28 dAbdo et al.ll2010bl). W44 dAbdo et al.l l20l0cT). IC 443 
dAcciari et dJ2009tlAbdo et al.ll2010dh G349.7+0.2. CTB 37A, 
3C391, and G8.7-0.1 dCastro & Slane1l2010l) are located close 
to dense regions or molecular clouds (MCs), thus one may use 
them to infer the diffusion parameters in thei r vicinity. Triggered 
by these o bservations, a number of studies dGabici et al.l l2009; 
iFuiita et aLll2009l) have attempted to explain the soft gamma-ray 
emission from MC-SNR systems in terms of escape d CRs. These 
studies have been based on ana lytical models dAtovan et al.l 
H99H lAh aronian & Atovani n~996) and several assumptions: (i) 
the SNR evolution is either stationary or Sedov-like, (ii) often 
that particle acceleration is quasi-instantaneous compared with 
the CR diffusion time, and (iii) the escaping-particle distribu- 
tion is a power-law or monoenergetic at a given maximum mo- 
mentum. The source of CRs is considered to be point-like, i.e., 
the SNR radius is much smaller than the distance to the MC . 
iLi & Chenl d2010h . iLi & Chenl d2011h . and lOhira et all (l2011h . 
extended the analytical models to include finite-size sources and 
finite-size target MCs. 

Despite their simplicity, the analytical studies demonstrate 
that to explain observations, the diffusion coefficient in the SNR 
vicinity must be roughly two orders of magnitude sma ller than 
the av erage Galactic value, supporting earlier claims dWentzell 
1974) that the diffusion coefficient might be smaller because 
of plasma waves that scatter particles. Using a simple model 
of SNR evolution and Monte Carlo simulations of CR diffu- 
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sion, Fuiita et al. (2010) showed that particles may be trapped 
around the SNR f or a significant period of time. Subsequently, 
Fuiita et al] d201 ll) confirmed this result using a transport equa- 
tion to describe particle acceleration at the forward shock (FS), 
and then a simplified transport equation to follow the propaga- 
tion of particles with p/moc > 10 2 4 from the precursor bound- 
ary (set at four Bohm diffusion length-scales) to distances far 
away from t he shock. These find ings are consistent with the con- 
clusions of iReville et alj d2009l) . who argue that beyond some 
critical value, L, which they call the free-escape boundary, the 
number density of CRs falls significantly. The position of L is 
dictated by the level of excited magnetohydrodynamics (MHD) 
turbule nce, which may be significant o ut to one SNR radius from 
the FS dZirakashvih & Pfuskinll2008h . The detailed variation of 
the turbulence amplitudes and the CR diffusion coefficient de- 
pends on the micro-phy sical balance of streaming instabilities 
and turbulence damping (lYan et al.ll2012l) . 

Besides the diffusion coefficient in the SNR vicinity, other 
important aspects of the study of escaped CRs are their spec- 
tral and spatial distributions. These distributions affect the ra- 
diative properties of the SNR itself as escaping CRs modify 
the high-energy part of the CR spectrum, and, if present, the 
radiative properties of the nearby dense matter consisting of a 
MC or a dense shell swept up by the winds of a high-mass SN 
progenitor). Much of the research dealing with stationary and 
(semi-)analytical solutions cannot recover the spectral and spa- 
tial shapes of the CR distribution, but only gives the integrated 
escaping CR energy flux. Kinetic models provide the shape of 
the distribution , which depends signific antly on the assumed dif- 
fusion models. Ellison & Bvkovl (1201 lb parametrized the spec- 
tral distribution of escaping CRs, and then propagated the es- 
caping CRs in the upstream region of the SNR obtaining spatial 
distributions. They used a spherically symmetric code for the 
hydro simulations of the SNR evolution that is identical to the 
code used in this work, but without the grid expansion that sig- 
nificantly increases the resolution. However, the particle accel- 
eration was treated in a plane -parallel steady-s t ate ap proxima- 
tion ignoring dilution effects. lEllison & Bvkovl d201 lb focused 
on the core-collapse explosion of a massive star in a low density 
medium surrounded by a dense shell, a case similar to our model 
for a core-collapse SNR (CC-SNR) in a wind bubble. 

In this paper, we investigate how both the spectrum of par- 
ticles escaped from the SNR into a nearby MC or shell and the 
particle emission due to interaction with dense matter depends 
on the acceleration history of young ejecta-dominated Type-la 
and Type-Ic (core-collapse) SNRs, using given different assump- 
tions about CR diffusion in the vicinity of the SNR. For this pur- 
pose, we combine a test-particle treatment of CR acceleration 
with an account of escaping particles. We solve the CR trans- 
port equation in a spherically symmetric geometry. Our calcula- 
tions are based on realistic high-resolution hydrodynamic sim- 
ulations of SNR evolution. We explore two SN types, focusing 
on the complexities in the circumstellar environment. In the case 
of Type-la SNRs, we assume that evolution proceeds in a uni- 
form ISM. We also consider a core-collapse SNR arising from a 
Wolf-Rayet (WR) progenitor, expanding in the wind-blown bub- 
ble created by the progenitor star. The SNR shock wave first 
evolves in the freely expanding wind of the progenitor, and then, 
beyond the wind-termination shock, in the shocked- wind region. 
Subsequently, we would expect the shock to impact the dense 
shell bordering the wind-blown bubble. However, our intention 
here is to study the interaction of the accelerated particles pro- 
duced by the SNR with this dense shell, hence we terminated 
the SNR evolution before the shock-shell impact took place. We 



considered particle acceleration at both the forward and reverse 
shock. In addition, in the CC-SNR case we consider particle re- 
acceleration at secondary shocks formed when the SNR forward 
shock collides with the termination shock of the wind zone. The 
spatial domain considered in our particle simulations extends to 
about 100 times the SNR radius to facilitate the study of parti- 
cle escape. The "absorbing" boundary located that far from the 
SNR does not affect the result of our simulations. It is the spatial 
profile of the diffusion coefficient in the vicinity of the SNR that 
determines the degree to which particles can escape from the CR 
precursor. To estimate the observational consequences of CR es- 
cape, we consider the interaction of the escaped particles with a 
target MC close to a Type-la SNR, or in the core-collapse case 
we consider the interaction of CRs with the dense shell. The 
MC/shell are located well inside the simulation domain, hence 
we know the cosmic-ray number density at their location at any 
time. Using two different models of the CR diffusion upstream 
of the FS, we calculated the evolutions of both the cosmic-ray 
spectra at the location of the MC/shell and the hadronic gamma- 
ray emission as the SNR shock approaches the MC/shell. If not 
stated otherwise, we refer to particles as escaped when they have 
escaped from the precursor region of the FS and are found at the 
location of the MC or shell. The respective particle spectra are 
averaged over the volume of the MC or shell. 

2. Method 

We describe the hydrodynamics of the expanding SNR by nu- 
merical simulations. For t he Type-la SNR, we use the e jecta- 
density profile described in Dwarkadas & Chevalierl dl998l) . ex- 
panding into a constant- density medium. The s imulation itself 
is described in detail in lTelezhinskv et ah! (|2012). The evolution 
of core-collapse Type-Ic SNRs is much more complicated. We 
assume that the f ormation of a wind-blow n bubble is similar to 
that described bv lDwarkadasI dlool l2007h . The WR wind is as- 
sumed to have a mass-loss rate of 10~ 5 M o yr~', and a velocity 
of 1500 km s . The wind termination shock is assumed to be a 
strong shock that forms at a radius of 7 pc. As the SNR shock 
expands outwards in radius, it evolves in the freely expanding 
wind, impacts the wind termination shock, and then evolves into 
a constant density shocked wind. Finally, although this is not in- 
cluded in the current simulation, the shock will collide with the 
dense wind-blown shell surrounding the bubble. 

We treat CRs as test particles in gas-flow profiles given by 
the si mulations described above. Our method (Telezh insky et al.1 
2012) is based on a numerical solution of the CR transport equa- 
tion in a grid co-moving with the shock wave. To ensure suffi- 
cient resolution near the shock, the spatial coordinate is substi- 
tuted with a new coordinate, x„, for which a uniform grid is used 
when solving the particle transport equation 

(x-l) = (^--lj = (x t -l) 3 . (1) 

Thus, a rather coarse grid in i, is transformed into a very fine 
grid in x close to the shock, where the high resolution is needed 
to properly account for the acceleration of newly injected par- 
ticles. At the same time, this transformation allows us to sig- 
nificantly extend the grid far into the ISM (x » 1) with only 
a small extension in x„. The resolution obviously deteriorates 
with distance from the shock, but the mean free path of the par- 
ticles in the CSM/ISM is orders of magnitude larger than in 
the shock vicinity, thus permitting the use of a moderate reso- 
lution. The extension of the grid to several dozens of SNR radii, 
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and the spherical geometry, imply that the volume upstream of 
the shock is much larger than downstream. The co-moving grid 
obliterates any need to introduce an artificial "absorbing" or "es- 
cape" boundary because the particles will never reach it during 
the simulation time. The number of particles escaping from the 
SNR shock is defined purely by the diffusion properties of the 
upstream region. We do not need to define the diffusive flux 
through some artificial boundary to count the particles escaping 
through this boundary. In the current approach, owing to the grid 
extension the total number of injected particles is conserved and 
distributed over the considered volume according to diffusion. 
By the end of the simulation, we observe at the boundary of the 
simulation domain only numerical noise, which indicates that 
the boundary is reached by an insignificant fraction of particles. 
Technically, a macroscopic treatment of diffusion does not dif- 
ferentiate between individual particles bouncing back and forth 
in the shock precursor, thus all particles found upstream of the 
FS can be assumed to have escaped from the SNR0. Here we are 
interested in particles found at the location of the MC/shell (suf- 
ficiently far away from the shock precursor), and so we refer to 
these particles as escaped. 

Since our method is based on the test-particle approxima- 
tion, it requires that the CR pressure at the shock be less than 
10% o f the ram pressu re. We use a thermal-leakage injection 
model (B lasi et alJ2 005). and adjust the injection so that the CR- 
pressure limit is not violated. The injection coefficient, a free pa- 
rameter, is adjusted to be approximately 3 • 10~ 7 for the Type-la 
SNR, and 5 ■ 1(T 6 for CC-SNR. 

3. Magnetic field and diffusion 

3.1. Magnetic field parametrization 

Observations of young SNe, for e xample SN 1993 J 

(IChandra et al.l I2004I) or Tycho's SNR (lAcciari et all l201lh . 
suggest that the magnetic field (MF) required to explain the 
observed level of emission is orders of magnitude stronger 
than the average ISM field. We therefore assume that the MF 
in the shock vicinity is amplified, e.g. by streaming CRs dBelll 
l200l . Although in general both resonant and non-resonant 
modes operate, non-resonant amplific ation probably dom inates 
in the early stages of SNR evolution (Capr ioli et al.l l2009). The 
amplified field is then given by 



(2) 



where p u {i) is the density upstream the shock, V s (t) is the shock 
speed, c is the speed of light, and £(f) is the ratio of cosmic- 
ray to ram pressure. In our calculations, £(?) « 0.05 with small 
deviations throughout the simulation. 

We parametrize the MF inside and outside the SNR. The 
scaling inside follows the time-dependent density distribution 



B(r, t) = cr 



B Q (t)p(r, t) 
p(R FS ,t) ' 



(3) 



where <x = VTT is the compression ratio of the turbulent MF, 
and is assumed to be isotropic. We assume that the MF falls off 
exponentially down to the strength of the interstellar field (5 fiG) 
(or circumstellar MF (CMF) in the case of core-collapse SNR), 



1 This is not true in a microscopic treatment, where any individual 
particle has a non-zero probability of returning back to the shock from 
any distance. 
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Fig. 1. Time dependence of magnetic-field profiles for type-la 
(top) and core-collapse (bottom) SNRs. 



at 0.05 ^ FS ahead of the FS (Zirakas hvili & Aharonianl l2010), 
and likewise down to the very small ejecta field (0.01-0. 1 fiG) at 
0.05 R RS toward the interior of the SNR. 

To define the magnetic field in the wind zone of the WR 
star, we assume that the dominant component at large distances 
is merely the toroidal component of the stellar surface magnetic 
field, which decreases outwards with radius. Thus, this compo- 
nent is defined as 



B C (R) 



B S R 



WR 



R 



(4) 



where B s w 100 G is the MF at the surface of WR star, R WR = 
8/? Q is the radius of WR star, and R is the distance from the star. 
Beyond the WR wind zone, in the constant density region, the 
MF is assumed to be constant and equal to the shock-compressed 
value at the edge of the WR wind zone. We assume a constant 
MF in the MC/shell, since it has been found that the MF does 
not increase if the n umber density in the cloud, n < 300cm 3 
(ICrutcher et al.ll20~Toh . The evolution of the MF profiles for two 
types of SNR is shown in Fig.Q] 

3.2. Diffusion models 

To study the propagation of cosmic rays in the vicinity of SNRs, 
it is necessary to make assumptions about the efficiency of dif- 
fusion close to the SNR shock. We note that the type of par- 
ticle scattering does not depend on the amplitude of the MF, 
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but rather on the wave-number spectrum of MHD turbulence. 
In other words, one can have Bohm-type diffusion in the non- 
amplified MF even significantly far away from the shock. Bohm 
diffusion is widely assumed to operate close to the shock and 
inside the SNR, otherwise DSA theories would fail to explain 
particles with energies up to the "knee". On the other hand, the 
average diffusion coefficient in the Galaxy is much larger than 
Bohmian, and inevitably there must be a transition between the 
small Bohmian coefficient and the large Galactic coefficient. 

Since the transition from Bohm diffusion in the vicinity 
of the forward shock to Galactic diffusion further away is not 
clearly understood, we consider two different models of CR dif- 
fusion outside the SNR. In model Dl, the diffusion coefficient is 
Bohmian inside the SNR and in the FS upstream region, up to 
a fiducial boundary, L - 2Rsnr, at which the transition to the 
Galactic diffusion occurs. We emphasize that the boundary L is 
not a physical boundary in any sense, and does not denote an 
"escape" boundary as in our previous work. It is merely a con- 
venient radius at which the transition in diffusion coefficient is 
assumed to occur. In this model the diffusion coefficient is de- 
fined as 



D(r) 



D B 



r<L, 
r > L, 



(5) 



D fi = 



(6) 



where the Bohm diffusion coefficient is given by 
pvc 
3qB 

while the Galactic d iffusion coefficient is taken to be 
dBerezinskii et alJI 1990b 



(7) 



where E is the CR energy and Do is the normalization. 
Follo wing a statistical analysis of cosmic -ray propagation mod- 
els dTrotta et al.ll201 lh . we use Do = 10 29 cm 2 /s and a - 1/3. 

In model D2, we explore a less abrupt transition to Galactic 
diffusion. We assume that Bohm diffusion is valid both inside 
the SNR and in the FS upstream region stretching out to 5% 
of the SNR radius, i.e., up to I - \.05Rsnr- Between I and L, 
we introduce an intermediate diffusion coefficient, smaller than 
the Galactic diffusion coefficient, thus mimicking the presence 
of MHD waves which may be invoked by CRs streaming away 
from the SNR dYan et alj|20l2h . We define the diffusion coeffi- 
cient in D2 model as 



D(r) 



D1 



C D B R SNR <r<l, 
\xD c l<r<L, 
[Dc r>L, 



(8) 
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Fig. 2. Schematic view of diffusion models. 



where x = 0.01. A schematic view of the diffusion models and 
corresponding boundaries is given in Fig. [2] 

We do not consider the confinement of particles by the 
dense shell at the boundar y of the WR win d zone in the CC- 
SNR models as done by (Ellison & Bvkovl d2011l) . The ion- 
neutral friction in dense and partially ionized media may sup- 
press p lasma instabilities, and thus increase the diffusion coef- 
ficient dPtuskin & Zirakashvilil l2005t lDrurvl|201 lh. Similar ar- 
gumen ts would apply to MCs. For instance. lEverett & Zw eibel 
(1201 lh consider the diffusion coefficient in MCs to be l%-500% 
of the Galactic value. They claim that this variation does not sig- 
nificantly affect their findings, which state that the cosmic -ray 
density does not increase within clouds. Therefore, we assume 
that the diffusion coefficient inside the MC/shell is the same as 
in the surrounding medium. A possible variation of the diffusion 
coefficient inside the MC/shell and the resulting consequences 
for our calculations are discussed in subsection l4.ll 

4. Results and discussion 

We have calculated the confined and escaped cosmic -ray spectra 
at 400, 1000, and 2000 years after explosion, assuming the diffu- 
sion models described above. We have also calculated the corre- 
sponding emission from the SNRs and the target material in their 
vicinity. In calculating the emission, we considered only CR pro- 
tons, and their radiation via pion decay (Hu ang et al.l l2007). 

Given the initial parameters assumed in the simulation, the 
Type-la SNR approaches the Sedov stage after around 1000 
years. Beyond this time, w e approximate the pl asma-flow pro- 
files with a Sedov solution (ICox & Francdfl98ll) . Although hy- 
drodynamically the RS is still present, usage of the Sedov so- 
lution is justified since the contribution from the RS to particle 
and e mission spectra is negligible at that time dTelezhinskv et al.l 
2012). After 2000 years, the FS radius in Type-la SNR is 
about 8 pc. The CC-SNR, on the other hand, by design, ex- 
pands in a wind of much lower density, followed by a constant- 
density shocked-wind zone, and therefore remains in the ejecta- 
dominated stage until the end of the simulation. At the edge of 
the constant-density zone, we assume the presence of a cold, thin 
(0.5 pc) dense shell of swe pt-up material, sim ilar to those known 
to exist around WR stars dCappa et alj|2003h . After 2000 years, 
the forward shock of this Type-Ic SNR has expanded out to a 
radius of about 14 pc. 

The spectrum of particles interacting with the MC is a func- 
tion of its distance from the SNR shock, as well as the diffusion 
coefficient. We consider two center-to-center distances from the 
Type-la SNR to the MC, 12 pc (the "near" scenario) and 16 pc 
(the "far" scenario), thus study the effect of a variation in the 
diffusion coefficient on the particle spectrum as the cloud ap- 
proaches the SNR. The MC is assumed to have a radius of 4 pc 
and a mass of M c m 1500 M , corresponding to a number den- 
sity n c =150 cirT 3 . 

For the CC-SNR case, we consider a swept-up shell located 
at 16 pc (the "near" scenario) or at 30 pc (the "far" scenario). 
The number density of the shell, n s = 100 cm 3 , corresponds to 
a mass M s * 1300 M for the "near" and M s ~ 4500 M for the 
"far" scenario. We note that the gamma-ray flux simply scales 
with the gas mass in the MC or the dense shell, and therefore the 
choice of mass only affects the normalization of the spectra. 



4. 1 . Particle spectra 

The time evolution of the spectra of confined and escaped pro- 
tons are presented in the upper panels of Fig.[3]for Type-la SNR 
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Fig. 3. Results for Type-la SNR. Top: the time evolution of volume-integrated CR spectra inside the SNR (denoted "snr") and the 
MC (denoted "mc") for two different diffusion models ("D1"/"D2") and two different distances from the SNR to MC: 12 pc ("n") 
and 16 pc ("f"). Bottom: the corresponding hadronic gamma-ray emission spectra. The thick solid lines represent the Galactic CR 
background and the corresponding gamma-ray emission from the MC, both of which do not depend on time. The thin solid lines 
represent the spectra of CRs outside the SNR for the D2 model and the corresponding emission from ambient diffuse gas. 



and Fig. [4] for CC-SNR, respectively. All spectra are volume- 
averaged. The Galactic CR background and its gamma-ray emis- 
sion are shown for comparison. 

We note from the figures that the choice of diffusion model 
influences the cut-off energy in the spectra of the confined parti- 
cles, as well as the energy of the spectral peak of escaped CRs. 
Particle acceleration becomes slow when the CR precursor ex- 
tends to the radius I in the upstream region at which diffusion 
becomes faster in the D2 model. Thus, the spectra of confined 
particles cut off at lower energies, and likewise the spectra of es- 
caped CRs have maxima at lower energies, because escape from 
the SNR is more efficient. In both models of diffusion and both 
types of SNR, we do not observe any significant contribution 
of the RS to the population of the CRs found upstream of the 
FS. This is partly because the reverse-shock contribution dimin- 
ishes with time, and is largely negligible after about 1000 years. 
Furthermore, since the MF is assumed to scale with the SNR 
density in the shocked interaction region, the large MF at the 
contact discontinuity, and the low maximum energy attained by 
CRs at the reverse shock do not permit particles accelerated at 
the RS to diffuse out of the SNR. 

In Type-la SNR, the spectral shapes of escaped particles in 
the Dl and D2 diffusion models are similar. They have a nar- 
row, parabola-like shape at a certain energy, usually referred 
to as ffmax, that i s consi stent with the approximations made by 
Ellison & Bvko^ (1201 ll) . The intensity changes with time de- 
pending on the assumed distance to the MC. In the "far" sce- 
nario, the MC always stays beyond the L boundary, and the in- 
tensity of the CR spectra at the location of the MC rapidly de- 
creases with time on account of dilution. In the "near" scenario, 
the MC is partially (after 1000 years) or totally (after 2000 years) 



within L, where the number density of CRs is high. However, the 
CRs in the D2 model have a shallower gradient than in the D 1 
model, and so the CR illumination of the MC increases even if 
the MC is only partially within L (1000 years). In the Dl model, 
the spectrum is still affected by dilution, and it is only later (af- 
ter 2000 years), when the MC is totally inside L, that the spectral 
intensity rises dramatically. 

The behavior of escaped particles in the Dl model is sim- 
ilar for both types of SNR. However, the spectral shape of the 
escaped particles in the D2 model is somewhat different in the 
CC-SNR case. The spectral peak for the escaped CRs shifts to 
higher energies with time. Therefore, an asymmetry is seen in 
their spectra, most visibly so at later times. 

We now discuss the consequences of a different diffusion co- 
efficient inside the MC/shell. First, we note that the shell thick- 
ness is so small that the shell spectra are virtually unaffected, un- 
less the diffusion coefficient inside is extremely small. Suppose 
then that the plasma waves are damped so effectively that the 
diffusion coefficient inside the MC increases significantly (about 
ten times) with respect to the diffusion coefficient of the sur- 
rounding medium. It is obvious that the CRs will quickly attain 
a uniform distribution inside the MC, and the CR number den- 
sity at the far side of the MC will be nearly equal to that at the 
near side. This should affect the CR intensity inside the MC af- 
ter about 1000 years, when the MC is located at a position cor- 
responding to a strong gradient in the CR distribution (it is par- 
tially within boundary L). Estimates show that the CR intensity 
in this case may increase by a factor of a few. At other times, the 
differences will be marginal. 

The diffusion coefficient inside the MC may also be lower 
than in the ISM. This may be due to a strong gradient in CRs at 
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Fig. 4. Results for core-collapse SNR. Top: the time evolution of volume-integrated CR spectra inside the SNR (denoted "snr") and 
at the dense shell (denoted "sh") for two different diffusion models ("D1"/"D2") and two different distances from the SNR to the 
shell: 16 pc ("n") and 30 pc ("f"). Bottom: the corresponding hadronic gamma-ray emission spectra. The thick solid lines represent 
the Galactic CR background and the corresponding gamma-ray emission from the dense shell, both of which still do not depend 
on time. The thin solid lines represent the spectra of CRs outside the SNR for the D2 model and the corresponding emission from 
ambient diffuse gas. 



the outskirts of the MC, which may trigger plasma instabilities 
and in turn reduce the CR diffusivity. A reasonable value could 
be O.OIDg- We note that this case may affect the Dl model after 
1000 years in the same manner as described above, because the 
diffusion coefficient inside the MC will be larger than within the 
boundary L. Spectra at other times and for other models should 
not be affected. The main effect of a reduced diffusion coefficient 
inside the MC is the longer penetration time and consequently a 
delay in the illumination of the MC. If one estimates the time 
needed for a 10-TeV particle to reach the center of the MC, 
one finds t - R 2 mc /2Dmc — 230 years, where Dmc — 0.0 IDg- 
However, the outer 1.5-pc thick shell of the MC would then sup- 
ply 80% of the flux, and a 10-TeV CR would need only -30 
years to penetrate to this depth. The delay in illumination is 
therefore negligible compared with the uncertainties in the size 
and structure of the MC. 



4.2. Spatial distributions of CRs 

To study the propagation of escaping particles in the medium 
around SNRs and understand the significance of their contribu- 
tion to the local CR level, we plot the radial distributions of par- 
ticles of given energies for different types of SNR and different 
diffusion model. Fig. [5] presents the radial distributions of pro- 
tons at the peak energy in the spectra of escaped CRs, which is 
sometimes referred to as E max . This energy changes with time, 
therefore at each time particles of different energy are plotted 
(thick lines), along with the corresponding Galactic background 
(thin lines). The evolution of £ mal is shown in Fig. [8] and dis- 



cussed below. We observe significant confinement of cosmic -ray 
particles close to the shock for all SNRs and diffusion models, 
even at the highest particle energies. Remarkably, for Type-la 
SNR in both diffusion models the highest-energy CRs still dom- 
inate over the background at a few SNR radii, whereas for CC- 
SNR this happens only in the D2 diffusion model. This may be 
because in our simulations the Type-la SNR propagates into a 
medium of higher average density, and therefore has a radius 
roughly twice as small as the CC-SNR, which implies less di- 
lution. The CC-SNR produce a lower CR energy density than 
Type-la SNRs owing to the lower density at the FS, hence lower 
injection into the acceleration process. Therefore, the intensity of 
CRs around CC-SNR only marginally exceeds the background. 
While these results are applicable to Type-Ic SNR, they are not 
to SNR of Type-IIP. The latter could have even smaller radii than 
Type-la, for instance, as they expand within the high-density red 
supergiant wind. In future papers, we will discuss the various 
SNR types in more detail. 

We are particularly interested in the spatial distribution of 
CRs that produce TeV-band gamma-ray emission. We show the 
time evolution of these distributions for a CR energy of 20 TeV 
in Fig. [6] including for comparison the Galactic CR background. 
We note from this figure that in the Dl diffusion model, 20-TeV 
CRs are trapped very close to the SNR shock (< IARsnr for 
Type-la SNR and < 1.2R SNR for CC-SNR) at all times. Even in 
the D2 model, in which the diffusion coefficient is very large, at 
late times the intensity of 20-TeV CRs exceeds the background 
at a few Rsnr on ly m the case of a Type-la SNR. Therefore, 
the detection of gamma-ray emission from MCs as probe of 
hadronic acceleration in SNR in the era of Fermi and current CTs 
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Fig. 5. Radial distributions of CRs at the peak energy in the spectra of escaped CRs at the given time (thick lines), compared with 
the Galactic CR background at the same energies (thin lines) for diffusion models "Dl" (left) and "D2" (right). Type-la SNRs are 
shown in the top row and core-collapse SNR in the bottom row. 



is a far more complicated task than previously thought. The ob- 
ject observed is always a composite of the SNR itself, emission 
from the MC, or a dense shell illuminated by freshly accelerated 
CRs, and the emission of interactions between upstream CRs 
and the ambient diluted gas. One may be successful in scenarios 
in which the MC is located close to the SNR. The distance to the 
MC should then be generally smaller or comparable to Rsnr, 
otherwise at the location of the MC/shell the background CRs 
dominate over the CRs from the SNR. This is a direct conse- 
quence of both diffusion and dilution, which becomes a serious 
factor at large distances from the SNR. 

4.3. Gamma-ray emission from the SNR-MC/shell 

We calculated the gamma-ray emission from the MC/shell il- 
luminated by CRs that had escaped from the SNR. For com- 
parison, we plot their spectra together with the emission spectra 
from the SNRs themselves (Fig.[3]and Fig.|U lower panels). The 
gamma-ray emission is only moderately affected by the choice 
of diffusion model. In the case of the D2 model, we observe 
higher fluxes at lower energies and also that the spectral peak 
is shifted to lower energies, thus making the spectra soft in the 
band above ~ 10 TeV. The spectra thus reflect the spatial vari- 
ation in the diffusion coefficient in the vicinity of the SNR. In 
Type-la SNR, the MC is always significantly brighter than ex- 
pected when illuminated by Galactic CRs. In contrast, the dense 
shell around a CC-SNR does not experience a significantly en- 
hanced CR illumination until the shock is very close. The emis- 
sion from the CC-SNR itself is never particularly bright, and 



so at late times the shell emission strongly dominates. On the 
other hand, the Type-la SNR dominates the intensity distribution 
of the SNR-MC system in all scenarios and diffusion models, 
but only at the latest time, and in the "near" scenario the emis- 
sion levels of SNR and MC become comparable. The intensity 
of gamma-ray emission from the MC/shell obviously strongly 
depends on the mass carried by the MC/shell, but the trend de- 
scribed above holds and may permit us to differentiate between 
Type-la and Type-Ic SNRs . 

Additionally, for model D2 we constructed the intensity dis- 
tributions of gamma-ray emission from Type-la SNR/MC and 
CC-SNR/shell systems, along with the background created by 
the upstream CRs in the vicinity of the SNR (see Fig. [7J. Since 
the number density of CRs confined within boundary L signif- 
icantly exceeds the number density of CRs beyond L, we con- 
sider L as the boundary of the SNR vicinity. One can see how 
the expanding halo of upstream CRs illuminates the surround- 
ing diluted gas and dense matter. The sharp contrast in bright- 
ness between the inner and outer parts of the MC at the age of 
1000 years arises from the abrupt transition from intermediate 
to Galactic diffusion and is in that sense an artifact. Moreover, 
as noted earlier in subsection 14. 11 if the diffusion coefficient in- 
side the MC were larger than that in the surrounding medium, 
the CR illumination of the MC after 1000 years would be more 
uniform so that the far side becomes brighter. However, the im- 
age, though idealized, is not unrealistic. The bright outer ring 
in CC-SNR intensity maps is the swept-up shell illuminated by 
escaping CRs. The smaller ring at the age of 1000 years corre- 
sponds to the contact discontinuity (CD) of the CC-SNR. The 
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Fig. 6. Radial distribution of CRs with energy 20 TeV at different times compared with the Galactic CR background at the same 
energy (dotted line) for diffusion models "Dl" (left) and "D2" (right). Type-la SNRs are shown in the top row and core-collapse in 
the bottom row. 



FS can almost not be seen owing to the very low density in the 
wind zone. At the age of 2000 years, the shell strongly domi- 
nates, and the shocks are not seen. There is only a marginal hint 
of emission from the CD region. 



Different authors may define E max in different ways, e.g., the 
maximum energy at which particles leave the system, the cut- 
off energy, or the peak in the spectrum of escaped particles. In 
practice, these definitions should provide nearly identical values. 
In some models of CR acceleration, including analytical studies 
of CR escape, E max is a free parameter used in lieu of a free- 
escape boundary. 

The time evolution of E max is a crucial issue. We derive E lmx 
from fitting the spectra of escaped CRs for each SNR type and 
diffusion mod els. The results are plo tted i n Fig. [8] along wit h 
the results of lEllison & Bvkovl d201 ll) and iGabici etall d2009h . 
Although we consider MF amplification (M FA), we do not fin d 
as rapid a decrease in E max as predicted by Gab ici et alJ (12 009). 
Our cu rves are more compatible with those of lEllison & Bvkovl 
d2011l) . although they did not include MFA. Most analytical 
studies are based on the Sedov solution, MFA, and the assump- 
tion that particles start leaving the system at the beginning of 
the Sedov stage w hen the highest energies are reached (e.g. 
Gabi ci et alJ l2009). When we accurately account for the evolu- 
tion of a young SNR, in which the shock velocity is not constant 
and falls with time, and we include the effects of dilution by con- 
sidering spherically symmetric systems, we find that particles 



leave the SNR during the "free expansion" stage. In addition, 
lDTurvl(l20Tll) reminds us that escape is a random process and it 
is only in a probabilistic sense controlled by the ratio of the dif- 
fusion length to SNR radius. Therefore, in our opinion solving 
one diffusion-advection equation for all CRs, both those trapped 
inside the SNR and those "escaping", is the best way to treat the 
problem. At any instant, one can easily see which particles are 
inside and which outside the SNR. Defining E max by the spec- 
trum of CRs that have escaped into the MC/shell and observing 
the time evolution of £ max , one notices a dramatic difference be- 
tween our result and that derived with Sedov scaling. Moreover, 
as shown in Fig. [8] besides the early SNR evolution and SNR 
type, the evolution of £ max may also depend on the assumed dif- 
fusion model. We must note, however, that the trend seen in CC- 
SNRs for a D2 diffusion model may be an artifact arising from 
our choice of transition boundary from Bohm to intermediate 
diffusion coefficient, x^c- in our case, the transition is at 5% of 
the SNR radius, which given the CC-SNR evolution, is roughly 
equal to the size of the CR precursor, DB(E mdx )/V s . If the CR 
precursor extends into a region with a large diffusion coefficient, 
shock acceleration loses efficiency, thus effectively prescribing 
the E max . With time and increasing SNR radius, the physical sep- 
aration of the forward shock and the transition point in the dif- 
fusion model increases, thus permitting particles to reach higher 
energies. 
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Fig. 7. 1-TeV Intensity distributions of type-Ia/MC (top) and CC-SNR/shell (bottom) at the age of 1000 (left) and 2000 (right) years. 
The log-scaled colormap spans roughly over 2.5 orders of magnitude for each image. 



4.5. General remarks 



Ellison & Bvko^ d201 ll) provided analytical approximations of 
the spectra of escaping particles, which were then incorporated 
into their Monte Carlo cosmic -ray propagation model. Our re- 
su lts on the spectra of es caped CRs are compatible with those 
of lEllison & Bvkovl j20TTh . if the shock is sufficiently far away 
from the location of interest and the system is still young. At 
later times, the time evolution of £ mal (see Fig. leads to the 
formation of tails in the spectra of escaped CRs, that are popu- 
lated by particles that left the remnant early (see the 2000-year 
plot at Fig. gj. 

We noted that for both types of SNR and diffusion model 
(especially in Dl) we find a significantly smaller intensity ratio 



of esc aped CRs to trapped CRs at E max than lEllison & Bvkovl 
(l20Tlh . even considering their use of NDSA. This discrepancy 
possibly arises from the neglect of dilution in their model. 

We compared the radial profiles with those presented in the 
literature. In general, the radial pr ofiles obtained in ou r D2 mod- 
els are compatible with those in iFuiita et all d201 ll) (model C 
di sregarding the growth of Alfven waves) and with the profiles 
of Ellison & B vkovl d201 ll) (propagation in pure CSM) out to a 
radius of /?snr from the shock. These authors, however, did not 
make any comparisons with the background level as their results 
are presented in arbitrary units. We were unable to find research 
that shows the profiles for pure Bohm diffusion out to a full SNR 
radius from the shock. 



9 



I. Telezhinsky et al.: Time-dependent escape of CRs from SNRs, and their Interaction with dense media 



1000 



„ 100 
> 

CD 

X 

to 
E 

LU 

10 




la, D1 
la, D2 
Ic, D1 
Ic, D2 
GAC 2009 
E&B 2011 



400 



Fig. 8. Time evolution of E„ 
of escaped CRs. 



1000 
t[yr] 



2000 



, the peak energy in the spectrum 



4.6. Observational implications 

The TeV-band intensity m ap of the young T ype-la SNR Tycho 
observed with VERITAS (lAcciari et al.ll201ll) shows off-center 
emission, which may correspond to a fluctuation. However, it 
may also indicate an interaction of escaping CRs with a nearby 
MC. If this is true, Fig. [3] suggests that in the GeV band emis- 
sion from the SNR itself should dominate. So, if the TeV-band 
emission indeed partly originates from an illuminated MC out- 
side the remnant, the component from the SNR itself must have 
a softer spectrum than suggested by a naive comparison of Fermi 
and VERITAS data. 

Another interesting example is the you ng TeV-bright CC - 
SNR RX J 17 1 3 .7-3 946 detected with Fermi (lAbdo et al.ll20TTh . 
Ellis on et al.1 (1201 lh applied their model to explain its high- 
energy emission via escaping CRs and found that this scenario is 
very unlikely. However, Fig.|4]indicates peculiar spectral transi- 
tions, in particular for the D2 model, that make multi-component 
spectra quite conceivable in which most of the GeV-band emis- 
sion arises from the SNR itself, whereas a substantial fraction of 
TeV-band emission originates from gas exterior to the remnant. 
More careful modeling is needed to test the applicability of this 
scenario to RX J1713.7-3946. 

Although our treatment extends only over the first 2000 years 
of SNR evolution, we note that the rather peculiar spectra of 
CC-SNRs found in our models have spectral breaks that appear 
similar to those observed from some ol der Fermi SNRs (e.g. 
ICastro & Slanell2010t lAbdo et al.ll2010dl) . In addition, the ob- 
served discrepancy in the locations of the TeV and (sub-)GeV 
emission in the SNR vicinity can be considered in line with 
our models. For instance, as seen from the lower panels of 
Fig. [3] and |4] the emission from the dense matter far away from 
the SNR may dominate at TeV energies and disappear at low 
(sub-)GeV energies. At the same time, the emission from the 
SNR itself would be stronger at GeV and weaker at TeV ener- 
gies. Altogether, this makes the emission "peak" move over the 
image with the changing energy of the observations. In some 
sense, a similar sce nario is put forward by the AGILE team 
(iGiuiiani et al.l 1201 lh for the regions around SNRs W28 and 
IC 443, to explain why in the GeV-band the peak-emission error 
box is separated from the TeV-emission error box observed by 
VERITAS and MAGIC. 



5. Conclusions 

We have modified ou r modeling technique developed in 
iTelezhinskv et al.l (120121) to examine the distribution of escaped 
CRs at distances far from SNRs. Our approach combines a re- 
alistic treatment of SNR evolution using hydrodynamic sim- 
ulations with test-particle calculations of CR acceleration, by 
solving the transport equation of cosmic rays in a spherically 
symmetric geometry. Both the reverse and forward shocks are 
included in these simulations. We also account for the re- 
acceleration at reflected shocks in core-collapse SNR. Our ap- 
proach inherently includes the effects of dilution due to the ex- 
pansion of the system, as well as the finite size of the SNR and 
any gas target illuminated by escaping CRs. It permits the re- 
construction of the spectra of escaped CRs at any given distance 
from the source. 

We have shown that the peak energy and intensity of escaped 
particles strongly depends on the efficiency of the diffusion in the 
vicinity of the SNR. If diffusion is Bohmian out to the boundary 
L — 2R snr, CRs are very efficiently confined to the SNR, and 
only the highest-energy particles are able to diffuse out to some 
distance from the SNR. Even though CRs at E max leave the SNR, 
they are still trapped at distances far less than one SNR radius 
from the FS. If the diffusion coefficient in the FS upstream re- 
gion is much larger than Bohmian, but not at the level of average 
Galactic diffusion, then CRs of lower energy are able to escape 
from the SNR. Consequently, one observes the broader spectra 
of escaped CRs with lower E m;lx and the cut off in the spectra of 
confined CRs is slower than exponential. 

It is possible but difficult to constrain the diffusion coefficient 
in the vicinity of the SNR using the emission from nearby dense 
material. The particles cannot generally propagate far from the 
shock. Thus, to study CR acceleration in SNRs via CR escape 
into a nearby MC or shell, one needs to find MCs or shells lo- 
cated very close, unless the diffusion coefficient in the immediate 
environment of the SNR is a significant fraction of the average 
Galactic one, in which case E m - dx is lower. We note that in all 
diffusion models, and especially in the "far" scenarios, one can 
clearly observe the striking effect of dilution, i.e. the decrease in 
the CR intensity upstream of the FS, introduced by the spher- 
ical geometry of the SNR. Ignoring the spherical geometry of 
the system is a serious mistake when one models the gamma- 
ray emission from MCs placed at significant distances from the 
SNRs. 

We find that Type-la SNR in our simulation has a higher CR 
energy density in its vicinity than Type-Ic core-collapse SNR. 
Therefore, Type-la SNR may have brighter pion-decay gamma- 
ray emission than Type-Ic SNR that is still evolving in the wind 
zones of its progenitor star. The emission of Type-la SNR is 
normally brighter than the MC emission until the SNR shock 
is very close to the MC. Only then is the gamma-ray flux from 
the MC comparable to that from the SNR. The emission from 
the CC-SNR and the wind-blown shell is initially near the back- 
ground level. At a later stage, the shell emission can significantly 
dominate over the SNR emission. Although the exact scenario 
depends on the MC/shell parameters, most importantly the tar- 
get mass, the described trend may serve to differentiate between 
these SNR types. 

Finally, we find that the commonly used approximation for 
the E m - dx evolution based on Sedov scaling is not reproduced 
when the evolution of the SNR shocks is computed accurately. 
Moreover, the £ max behavior critically depends on the diffusion 
coefficient in the vicinity of the SNR. 
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